Read and organize data

read_srm_export <- function(filename, columns = c("peak_name", "RT.min", "basepeak", "area.cpm", "height.cts", "quantitation")) {
  filename %>% 
    # read excel files
    read_excel(sheet = "Integration", skip = 42, 
               col_names = columns, col_types = rep("text", length(columns))) %>% 
    as_data_frame() %>%
    # remove empty rows
    filter(!is.na(peak_name), peak_name != "n.a.") %>% 
    # convert the relevant numeric columns into numbers
    mutate_at(vars(RT.min, area.cpm, height.cts), as.numeric) %>% 
    # remove useless columns
    select(-basepeak, -quantitation) %>% 
    # add filename info
    mutate(file_id = gsub("\\.xls", "", basename(filename))) %>% 
    select(file_id, everything())
}

# get data
all_data <- 
  # find all excel files
  list.files("data_reports", recursive = TRUE, full.names = TRUE, pattern = "\\.xls$") %>% 
  # send them to the read method
  lapply(read_srm_export) %>% 
  # combine the data set
  bind_rows() %>% 
  # pull out sample information
  #mutate(sample_id = str_match(all_data$file_id, "TSQ\\d+_GB_(.*)$") %>% { .[,2] }) %>% 
  # get n replicates
  group_by(file_id)
  #mutate(n_replicates = length(unique(file_id)))

File names for metadata file

Calculation peak amounts and rock concentrations

depth_and_rock_info <- read_excel(file.path("metadata", "aromaticSRM_20180327.xlsx")) %>% 
  rename(rock.g = `rock .g`) %>% 
  filter(!is.na(file_id)) 
kable(depth_and_rock_info)
file_id depth rock.g process TLE.mg maltenes.mg
TSQ1998_GB_OG043 122.900 11.349 yes 14.7 11.9
TSQ2001_GB_OG058 121.195 10.711 yes 18.4 4.2
TSQ2002_GB_OG063 122.390 10.116 yes 139.5 8.8
TSQ2003_GB_OG043 122.900 11.349 yes 14.7 11.9
TSQ2004_GB_OG053 122.805 10.786 yes 71.2 6.7
TSQ2005_GB_OG054 122.700 12.054 yes 77.3 13.5
TSQ2006_GB_OG050 120.600 11.195 yes 16.1 2.1
TSQ2008_GB_OG048 121.295 10.828 yes 14.5 4.4
TSQ2009_GB_OG062 120.950 10.341 yes 83.1 7.8
TSQ2010_GB_OG011 114.000 10.208 yes 0.4 6.9
TSQ2011_GB_OG012 113.000 10.552 yes 9.4 11.4
TSQ2012_GB_OG013 111.080 9.747 yes 6.9 9.5
TSQ2013_GB_OG037 121.600 10.530 yes 5.2 3.6
TSQ2014_GB_OG038 120.700 11.444 yes 3.3 2.2
TSQ2015_GB_OG008 118.000 9.698 yes 3.4 5.0
TSQ2016_GB_OG040 121.890 10.300 yes 6.5 3.3
TSQ2018_GB_OG041 120.405 11.393 yes 2.3 1.9
TSQ2019_GB_OG042 120.205 10.757 yes NA NA
TSQ2020_GB_OG045 120.205 10.757 yes 4.1 2.4
TSQ2021_GB_OG044 120.300 11.390 yes 3.5 2.7
TSQ2022_GB_OG046 121.500 11.191 yes 24.7 3.1
TSQ2023_GB_OG047 121.105 10.868 yes 38.8 5.5
TSQ2024_GB_OG051 120.900 14.166 yes 96.8 6.3
TSQ2025_GB_OG052 121.800 10.268 yes 4.5 3.0
TSQ2026_GB_OG057 122.500 10.939 yes 36.8 12.6
TSQ2028_GB_OG059 122.090 9.486 yes 54.9 3.6
TSQ2029_GB_OG060 122.200 11.367 yes 53.3 2.8
TSQ2030_GB_OG061 121.400 12.044 yes 40.7 4.0
TSQ2031_GB_OG015 109.100 10.984 yes 13.3 15.3
TSQ2032_GB_OG007 119.000 9.620 yes -214.2 4.8
TSQ2033_GB_OG010 115.000 11.070 yes 5.0 6.0
TSQ2034_GB_OG009 117.010 10.742 yes 5.0 6.3
TSQ2035_GB_OG003 124.000 10.120 yes 224.4 7.3
TSQ2036_GB_OG005 122.000 9.280 yes 3.2 5.4
TSQ2038_GB_OG002 125.015 10.001 yes 24.9 13.7
TSQ2039_GB_OG039 120.100 7.997 yes 0.8 0.9
TSQ2044_GB_OG016 108.300 10.089 no 5.1 NA
TSQ2045_GB_OG017 106.900 10.348 yes -6.2 6.8
TSQ2046_GB_OG018 106.090 10.337 yes 4.9 6.5
TSQ2047_GB_OG019 105.100 11.002 yes 2.2 -5.2
TSQ2048_GB_OG021 103.890 10.056 yes 21.6 7.4
TSQ2049_GB_OG023 101.100 11.294 yes 8.5 NA
TSQ2051_GB_OG024 100.090 9.879 yes 6.3 5.6
TSQ2052_GB_OG025 98.900 10.713 yes 7.8 8.6
TSQ2053_GB_OG026 98.090 9.336 yes 4.7 7.4
TSQ2054_GB_OG027 96.890 9.378 yes 5.3 6.3
TSQ2055_GB_OG028 96.090 9.758 yes 6.9 7.3
TSQ2056_GB_OG029 95.080 11.469 yes 4.3 2.3
TSQ2057_GB_OG030 94.560 10.195 yes 7.1 8.2
TSQ2058_GB_OG031 119.950 7.502 yes 2.6 4.7
TSQ2060_GB_OG033 112.070 7.287 yes 2.3 -7.2
TSQ2061_GB_OG022 101.900 11.025 yes 4.5 7.2
TSQ2062_GB_OG014 110.100 10.273 yes 14.5 14.0
TSQ2063_GB_OG032 115.885 9.032 yes -2.1 4.4
137 117.800 11.205 yes 14.0 16.4
138 123.600 10.512 yes 15.5 31.6
139 118.925 10.663 yes 8.3 8.6
140 124.600 11.468 yes 12.3 7.6
141 118.490 10.656 yes 8.1 7.3
TSQ2221_GB_OG143 119.100 10.763 yes 284.8 5.0
TSQ2222_GB_OG144 116.700 10.862 yes 10.6 9.5
TSQ2223_GB_OG145 117.090 11.360 yes 9.5 2.8
146 NA 0.000 no 4.0 0.3
TSQ2224_GB_OG147 SB NA 0.000 no 3.6 0.3
TSQ2225_GB_OG148 117.210 10.792 yes 8.4 3.2
TSQ2226_GB_OG149 117.700 10.463 yes 6.3 2.0
TSQ2228_GB_OG150 117.700 10.393 yes 5.5 1.7
TSQ2229_GB_OG151 117.900 11.828 yes 5.1 1.0
TSQ2230_GB_OG152 117.280 11.463 yes 10.6 3.9
TSQ2231_GB_OG153 116.810 11.177 yes 9.8 3.9
TSQ2232_GB_OG154 116.600 12.080 yes 29.1 3.6
TSQ2233_GB_OG155 119.850 11.042 yes 7.1 1.6
TSQ2235_GB_OG156 117.400 11.244 yes 8.0 3.3
TSQ2236_GB_OG157 118.800 11.346 yes 8.3 2.6
TSQ2237_GB_OG158 119.550 11.853 yes 7.6 3.2
TSQ2238_GB_OG159 119.655 11.044 yes 6.7 3.8
TSQ2240_GB_OG160 118.300 11.305 no 7.4 2.1
TSQ2242_GB_OG161 116.490 10.770 yes 8.4 2.9
TSQ2243_GB_OG162 116.900 10.931 no 10.9 4.1
TSQ2244_GB_OG163 119.445 11.636 no 1.8 1.4
TSQ2215_GB_OG164 118.100 10.619 yes 7.6 3.5
TSQ2216_GB_OG165 119.570 11.300 yes 5.0 1.6
TSQ2217_GB_OG166 118.600 10.790 yes 7.1 1.9
TSQ2218_GB_OG167 117.490 10.853 yes 8.8 2.9
TSQ2219_GB_OG168 118.965 11.259 yes 6.3 2.0
data_by_depth <- 
  all_data %>%
  left_join(depth_and_rock_info, by = "file_id") %>% 
  group_by(file_id) %>% 
  mutate(
    n_peaks = n(),
    n_standards = sum(peak_name == "d14-pTerph"),
    ref_area.cpm = area.cpm[peak_name == "d14-pTerph"],
    ref_amount_added.ug = 5, #check standard amount added
    amount.ug = area.cpm/ref_area.cpm * ref_amount_added.ug,
    conc_rock.ug_g = amount.ug / rock.g
    #total_area.cpm = sum(area.cpm[peak_name != "d14-pTerph"]),
    #area.percent = area.cpm / total_area.cpm * 100
  )%>% ungroup() %>% 
  filter(`process` == "yes") %>%
  arrange(file_id, peak_name)

data_by_depth
## # A tibble: 4,705 x 16
##    file_id    peak_name    RT.min area.cpm height.cts depth rock.g process
##    <chr>      <chr>         <dbl>    <dbl>      <dbl> <dbl>  <dbl> <chr>  
##  1 TSQ2001_G… 1-MP           26.1   9.69e⁶  168777000   121   10.7 yes    
##  2 TSQ2001_G… 2-MP           26.8   1.26e⁷  228207000   121   10.7 yes    
##  3 TSQ2001_G… 3-MP           26.7   1.47e⁷  265019000   121   10.7 yes    
##  4 TSQ2001_G… 9-MP           25.9   1.09e⁷  198920000   121   10.7 yes    
##  5 TSQ2001_G… Acenapthene    14.2   2.71e⁶   56677300   121   10.7 yes    
##  6 TSQ2001_G… Benzo[a]ant…   41.8   6.84e⁶  105420000   121   10.7 yes    
##  7 TSQ2001_G… Benzo[a]pyr…   51.8   1.51e⁷  243142000   121   10.7 yes    
##  8 TSQ2001_G… Benzo[b]flo…   49.7   1.25e⁷  146978000   121   10.7 yes    
##  9 TSQ2001_G… Benzo[e]pyr…   51.4   1.84e⁷  300891000   121   10.7 yes    
## 10 TSQ2001_G… Benzo[ghi]     60.0   1.66e⁸ 1312680000   121   10.7 yes    
## # ... with 4,695 more rows, and 8 more variables: TLE.mg <dbl>,
## #   maltenes.mg <dbl>, n_peaks <int>, n_standards <int>,
## #   ref_area.cpm <dbl>, ref_amount_added.ug <dbl>, amount.ug <dbl>,
## #   conc_rock.ug_g <dbl>

Calculate Recovery

Linear regressions of the calibration curves

standard <- read_excel(file.path("metadata", "D14 calibration.xlsx"))   ###read excel

###calibration curve
standard %>% 
  ggplot() +
  aes(x = Known.ng, y = Measured_area.counts, color = calibration) + 
  geom_smooth(method = "lm", alpha = 0.5) +
  geom_point() +
  theme_bw() +
  theme(legend.position = "none") 

calibrations <- 
  standard %>% 
  filter(!is.na(calibration)) %>% 
  nest(-calibration) %>% 
  mutate(
    fit = map(data, ~summary(lm(`Measured_area.counts`~ `Known.ng`, data = .x))),
    coefficients = map(fit, "coefficients"),
    intercept = map_dbl(coefficients, `[`, 1, 1),
    intercept_se = map_dbl(coefficients, `[`, 1, 2),
    slope = map_dbl(coefficients, `[`, 2, 1),
    slope_se = map_dbl(coefficients, `[`, 2, 2),
    r2 = map_dbl(fit, "r.squared")
  )

calibrations %>% select(-data, -fit, -coefficients) %>% knitr::kable(d = 3)
calibration intercept intercept_se slope slope_se r2
sept2017 -131054.9 2131453 2275077 90941.37 0.997
oct2017 -1488460.7 2799385 1557777 107323.96 0.995

Determine yield

These numbers are not useful for anything else.

calib_data <-
  data_by_depth %>% 
  # temp
  mutate(calibration = "oct2017") %>% 
  left_join(calibrations, by = "calibration") %>% 
  mutate(
    total_volume.uL = 100,
    total_inject.uL = 1,
    ref_amount_inject_expected.ng = ref_amount_added.ug/total_volume.uL * total_inject.uL * 1000,
    ref_amount_inject_measured.ng = (ref_area.cpm - intercept)/slope,
    ref_amount_measured.ug = total_volume.uL/total_inject.uL * ref_amount_inject_measured.ng * 1/1000,
    yield = ref_amount_inject_measured.ng/ref_amount_inject_expected.ng
  )

check yields

calib_data %>% 
  select(file_id, yield)  %>% 
  arrange(file_id)  %>% 
  unique() %>% 
  ggplot() + aes(file_id, y = 100*yield) +
  geom_point(size = 3) +
  theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5))

Combine compounds/make ratios (new rows w/o RTs etc, just concentration.rock column)

# functions to make it easy to sum up peaks

sum_peaks <- function(df, filter_condition, new_peak_name) {
  filter_condition <- sprintf("(%s)", str_c(filter_condition, collapse = "|"))
  filter(df, str_detect(peak_name, filter_condition)) %>% 
    summarize(
      file_id = file_id[1],
      depth = depth[1],
      conc_rock.ug_g = sum(conc_rock.ug_g)
    ) %>% 
    mutate(peak_name = new_peak_name)
}

ratio_peaks <- function(df, filter_top, filter_bottom, new_peak_name) {
  filter_top <- sprintf("(%s)", str_c(filter_top, collapse = "|"))
  filter_bottom <- sprintf("(%s)", str_c(filter_bottom, collapse = "|"))
  filter(df, str_detect(peak_name, filter_top) | str_detect(peak_name, filter_bottom)) %>% 
    summarize(
      file_id = file_id[1],
      depth = depth[1],
      ratio = sum(conc_rock.ug_g[str_detect(peak_name, filter_top)]) / sum(conc_rock.ug_g[str_detect(peak_name, filter_bottom)])
    ) %>% 
    mutate(peak_name = new_peak_name)
}
final_data <- data_by_depth %>% 
    group_by(file_id) %>% 
        do({
          bind_rows(., 
              sum_peaks(., "C15", "all C15"),
              sum_peaks(., "C16", "all C16"),
              sum_peaks(., "C17", "all C17"),
              sum_peaks(., "C18", "all C18"),
              sum_peaks(., "C19", "all C19"),
              sum_peaks(., "C20", "all C20"),
              sum_peaks(., "C21", "all C21"),
              sum_peaks(., "C24", "all C24"),
              sum_peaks(., "C26", "all C26"),
              sum_peaks(., "Aryl", "all_Aryl_Isop"),
              #sum_peaks(., "MP", "3ring_MP"),
              sum_peaks(., c("Acenapthene", "Flourene"), "2ring_all" ),
              #sum_peaks(., "Acenapthene", "2ring_1"),
              #sum_peaks(., "Flourene", "2ring_2"),
              sum_peaks(., c("Phenanthrene", "Flourantherene", "Retene", "1-MP", "2-MP", "3-MP", "9-MP"), "3ring_all" ),
              #sum_peaks(., "Phenanthrene", "3ring_1"),
              #sum_peaks(., "Flourantherene", "3ring_2"),
              #sum_peaks(., "Retene", "3ring_3"),
              sum_peaks(., c("Pyrene", "Benzo[a]anthracene", "Triphenylene", "Chrysene", "flouranthrene"), "4ring_all" ),
              #sum_peaks(., "Pyrene", "4ring_1"),
              #sum_peaks(., "Benzo[a]anthracene", "4ring_2"),
              #sum_peaks(., "Triphenylene", "4ring_3"),
              #sum_peaks(., "Chrysene", "4ring_4"),
              #sum_peaks(., "flouranthrene", "4ring_5s"),
              sum_peaks(., c("Benzo[e]pyrene", "Benzo[a]pyrene", "Perylene", "Ideno[c,e]", "Dibenz[a,h]"), "5ring_all" ),
              #sum_peaks(., "Benzo[e]pyrene", "5ring_1"),
              #sum_peaks(., "Benzo[a]pyrene", "5ring_2"),
              #sum_peaks(., "Perylene", "5ring_3"),
              #sum_peaks(., "Ideno[c,e]", "5ring_4"),
              #sum_peaks(., "Dibenz[a,h]", "5ring_5"),
              sum_peaks(., c("Benzo[ghi]", "Coronene"), "6ring_all" ), 
              sum_peaks(., c("Acenapthene", "Flourene", "Phenanthrene", "Flourantherene", "Retene", "Pyrene", "Benzo[a]anthracene", "Triphenylene", "Chrysene", "flouranthrene", "Benzo[e]pyrene", "Benzo[a]pyrene", "Perylene", "Ideno[c,e]", "Dibenz[a,h]", "Benzo[ghi]", "Coronene"), "all_PAH")
              #sum_peaks(., "Benzo[ghi]", "6ring_1"),
              #sum_peaks(., "Coronene", "6ring_2")
 ) }) %>% ungroup() 
  
final_data
## # A tibble: 5,873 x 16
##    file_id    peak_name    RT.min area.cpm height.cts depth rock.g process
##    <chr>      <chr>         <dbl>    <dbl>      <dbl> <dbl>  <dbl> <chr>  
##  1 TSQ2001_G… 1-MP           26.1   9.69e⁶  168777000   121   10.7 yes    
##  2 TSQ2001_G… 2-MP           26.8   1.26e⁷  228207000   121   10.7 yes    
##  3 TSQ2001_G… 3-MP           26.7   1.47e⁷  265019000   121   10.7 yes    
##  4 TSQ2001_G… 9-MP           25.9   1.09e⁷  198920000   121   10.7 yes    
##  5 TSQ2001_G… Acenapthene    14.2   2.71e⁶   56677300   121   10.7 yes    
##  6 TSQ2001_G… Benzo[a]ant…   41.8   6.84e⁶  105420000   121   10.7 yes    
##  7 TSQ2001_G… Benzo[a]pyr…   51.8   1.51e⁷  243142000   121   10.7 yes    
##  8 TSQ2001_G… Benzo[b]flo…   49.7   1.25e⁷  146978000   121   10.7 yes    
##  9 TSQ2001_G… Benzo[e]pyr…   51.4   1.84e⁷  300891000   121   10.7 yes    
## 10 TSQ2001_G… Benzo[ghi]     60.0   1.66e⁸ 1312680000   121   10.7 yes    
## # ... with 5,863 more rows, and 8 more variables: TLE.mg <dbl>,
## #   maltenes.mg <dbl>, n_peaks <int>, n_standards <int>,
## #   ref_area.cpm <dbl>, ref_amount_added.ug <dbl>, amount.ug <dbl>,
## #   conc_rock.ug_g <dbl>

Bring in other data for plots

osisotope <- read_excel(file.path("metadata", "SH1_Osi_forGarrett.xlsx")) %>% 
  rename(depth = `Depth (m)`) 
# Bring in SH1 data from Jones et al., (2018)
cisotope <- read_excel(file.path("metadata", "Appendix_Table1_geochemistry.xlsx")) %>% 
  #rename columns
  rename(depth = `Abs. depth (m)` , d13c_org = `Average δ13Corg (‰ VPDB)` , carb = `%Carbonate`, TOC = `%TOC`, d13c_carb = `Average δ13Ccarb (‰ VPDB)`) %>%
  #remove columns not of interest
  select(-`stdev δ13Corg`, -`δ13Ccarb stdev`, -`d18O-avg`, -`d18O stdev`, -`∆13C`, -`d13c_carb`)
carb <- cisotope %>%
  ggplot() +
  aes(depth, carb) +
  geom_point() +
  geom_line() +
  scale_x_reverse() +
  coord_flip() +
  ggtitle("carbonate%")


TOC <- cisotope %>%
  ggplot() +
  aes(depth, TOC) +
  geom_point() +
  geom_line() +
  scale_x_reverse() +
  coord_flip() +
  ggtitle("TOC%") +
  theme(axis.text.y = element_blank(), 
            axis.ticks.y = element_blank(), 
            axis.title.y = element_blank())


d13c <- cisotope %>%
  ggplot() +
  aes(depth, d13c_org) +
  geom_point() +
  geom_line() +
  scale_x_reverse() +
  coord_flip() +
  ggtitle("d13C_org")

os <- osisotope %>%
  ggplot() +
  aes(x = depth, y = Osi) +
  geom_point() +
  geom_line() +
  #facet_wrap(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip() +
  ggtitle("Os") +
  theme(axis.text.y = element_blank(), 
            axis.ticks.y = element_blank(), 
            axis.title.y = element_blank())

#ggplotly(carb) 
#ggplotly(TOC)
#ggplotly(d13c)

library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
litho <- grid.arrange(carb, TOC, d13c, os, ncol=4)
## Warning: Removed 3 rows containing missing values (geom_point).

By compound

Carotenoids against depth

PZE <- subset(final_data, peak_name %in% c('Chlorobactane', 'Isorenieretane', 'all_Aryl_Isop')) %>%
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_wrap(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip() +
  ggtitle("PZE") 

ggplotly(PZE)
#litho <- grid.arrange(d13c, os, PZE, ncol=3)

all PAH combos

all <- subset(final_data, peak_name %in% c('all_PAH', '2ring_all', '3ring_all', '4ring_all', '5ring_all', '6ring_all' )) %>%
  #filter(!is.na(depth==108)) %>%
    ggplot() +
    aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
    geom_point() +
    facet_grid(~peak_name, scales = "free") +
    scale_x_reverse() +
    coord_flip()

ggplotly(all)

Individuals

twos <- subset(final_data, peak_name %in% c("Acenapthene", "Flourene")) %>%  
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_wrap(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip()

ggplotly(twos)
trees <- subset(final_data, peak_name %in% c("Phenanthrene", "Flourantherene", "Retene", "1-MP", "2-MP", "3-MP", "9-MP")) %>%  
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_grid(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip()

ggplotly(trees)

4 ring PAHs

four <- subset(final_data, peak_name %in% c("Pyrene", "Benzo[a]anthracene", "Triphenylene", "Chrysene", "flouranthrene")) %>%
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_grid(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip()

ggplotly(four)

5 ring PAHs

five <- subset(final_data, peak_name %in% c("Benzo[e]pyrene", "Benzo[a]pyrene", "Perylene", "Ideno[c,e]", "Dibenz[a,h]")) %>%
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_grid(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip()

ggplotly(five)

6 ring PAHs

six <- subset(final_data, peak_name %in% c("Benzo[ghi]", "Coronene")) %>%
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_wrap(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip()

ggplotly(six)

all individual PAHs

allind <- subset(final_data, peak_name %in% c("Acenapthene", "Flourene", "Phenanthrene", "Flourantherene", "Retene", "1-MP", "2-MP", "3-MP", "9-MP", "Pyrene", "Benzo[a]anthracene", "Triphenylene", "Chrysene", "flouranthrene", "Benzo[e]pyrene", "Benzo[a]pyrene", "Perylene", "Ideno[c,e]", "Dibenz[a,h]", "Benzo[ghi]", "Coronene")) %>%
  ggplot() +
  aes(x = depth, y = conc_rock.ug_g, color = peak_name) +
  geom_point() +
  geom_line() +
  facet_grid(~peak_name, scales = "free") +
  scale_x_reverse() +
  coord_flip()

ggplotly(allind)